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ABSTRACT 

Context. Low-mass X-ray binaries hosting neutron stars (NS) exhibit thermonuclear (type-I) X-ray bursts, which are powered by 
unstable nuclear burning of helium and/or hydrogen into heavier elements deep in the NS “ocean”. In some cases the burning ashes 
may rise from the burning depths up to the NS photosphere by convection, leading to the appearance of the metal absorption edges in 
the spectra, which then force the emergent X-ray burst spectra to shift toward lower energies. 

Aims. These effects may have a substantial impact on the color correction factor f c and the dilution factor w, the parameters of the 
diluted blackbody model F E at wB E (f c T cS ) that is commonly used to describe the emergent spectra from NSs. The aim of this paper 
is to quantify how much the metal enrichment can change these factors. 

Methods. We have developed a new NS atmosphere modeling code, which has a few important improvements compared to our 
previous code required by inclusion of the metals. The opacities and the internal partition functions (used in the ionization fraction 
calculations) are now taken into account for all atomic species. In addition, the code is now parallelized to counter the increased 
computational load. 

Results. We compute a detailed grid of atmosphere models with different exotic chemical compositions that mimic the presence of 
the burning ashes. From the emerging model spectra we compute the color correction factors f c and the dilution factors w that can 
then be compared to the observations. We find that the metals may change f. by up to about 40%, which is enough to explain the 
scatter seen in the blackbody radius measurements. 

Conclusions. The presented models open up the possibility for determining NS mass and radii more accurately, and may also act as 
a tool to probe the nuclear burning mechanisms of X-ray bursts. 
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1. Introduction 

Low-mass X-ray binaries (LMXB) that consist of a neutron 
star (NS) primary, and a secondary star less massive than the 
Sun, may exhibit thermonuclear ( type-I) X-ray bursts (see e.g. 
Lewin et al. 119931 : fStrohmaver & BildstenlbOOfl for review). X- 
ray bursts are produced when the mass accretion rate onto 
the NS is in a regime where unstable nuclear fusion burns 
the accumulated hydrogen and/or helium into heavier elements 
dWooslev & Taamlll976c iFuiimoto et al.lll981l) . The heat gener¬ 
ated in the burning js released as X-ray radiation from the NS 
photosphere (tJossI 119781) . and some X-ray bursts are so ener¬ 
getic that the luminosity exceeds the local Eddington limit. In 
these cases the radiation pressure lifts the photosphere from 
the NS surface, and we can obse rve so-c alled photospheric ra¬ 
dius expansion (PRE) bursts dHoffman et al.|[T980l) . These PRE 
bursts are particularly interesting as they can be used to con¬ 
strain the NS mass and radius by comparing the observed cooling 
tracks to the predictions from theoretical NS atmosph ere models 
dSuleimanov et ai]l201 1 hi 120 1 2L [Poutanen et alJl2014l) . 

The nuclear burning region is a thin layer l ocate d 
deep in the neutron star “ocean” dFuiimoto et alJ 1 198 ll: 
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iFushiki & Lambl 19871). a t column de nsity on the order of 
10 8 g cm -2 (iCumming & Bildstenli2000l ). Detailed one-zone nu- 
cleosynthetic studies and hydrodynamical models coupled with 
vast reaction networks show that large quantities of heavier ele- 
ment s are p roduced during the nuclear burning process (see e.g. 
iParikh et al.120131 for a review). In PRE bursts the conditions are 
favorable for large scale conv ectio n and the burning ashes may 
rise up to the photosphere dWeinberg et al.1120061) . The presence 
of these heavier elements in the H and/or He plasma can then 
significantly alter the properties of the atmosphere. 


Theoretical calculations show that the emergent energy spec¬ 
trum emanating from the hot X-ray bursting NS photosphere is 
closely approximated by a diluted blackbody Fe ~ u!B E (f c T c {\), 
where w is the dilution factor, B E is the Planck function, f Q = 
7/./7j.fr is the color correction factor, T e g is the effe ctive t emp er¬ 
ature, and T c is the observed color temperature dLondon et al.l 
[1984^ Suleimanov et al.l201 lbi. 12012) . If the burning ashes reach 
the NS photosphere, strong absorption edges may appear on 
top of the diluted blackbody spectrum. The strength of the 
edges depends on the composition of the ashes , bu t also on the 
NS photospheric t emp erature ( Weinberg et al.1120061) . Encourag¬ 
ingly, lin’t Zand & Weinbergl ( 20101) have found evidence that 
several spectra of PRE bursts with very strong superexpansion 
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show strong residuals that can be modeled as absorption edges 
of heavy elements like 56 Ni. However, this evidence is not con¬ 
clusive with the current instruments, which have poor energy 
resolution. The origin of the absorption edges are hard to ver¬ 
ify, for example, the edges (and the emission lines) may also be 
produced b^ reflection from the surrounding accretion disk (e.g. 
IStrohmaver & Brownll2002h . 

Fortunately, the burning ashes have another effect on the 
X-ray burst spectra. If the metals reach the photosphere, they 
will increase the bound-free and the free-free opacity, and also 
the photospheric electron density, causing the emergent spec¬ 
tra to shift toward lower energies. In this case, the color cor¬ 
rection factor f c may decrease substantially. Given that each X- 
ray burst is likely to ignite in slightly different initial conditions, 
the amount and the composition of the burning ash that reaches 
the photosphere is expected to be variable between different 
bursts. Therefore, the burning ashes may also manifest them¬ 
selves as small burst-to-burst variations in the cooling tracks of 
individual bursts because the observed blackbody radii have a 
strong dependency on the color correction factor (/?bb x ff 2 )- 
In LMXBs where the RXTEfPCA have detected several PRE 


bursts from one source, one indeed sees such Rhh-scatter (e.g . 
Bhattacharvva et ai]l2010l : lzhang et al1l201 ltlGuver et al.ll2012t 


Poutanen et al .1120141) . Even if there are other possible mecha¬ 

nisms that may be responsible for it (such as anisotropic burst 
emission), the fluctuations of f c by the burning ashes is one of 
the strongest candidates. Motivated by this possibility, we have 
undertaken a study to quantify the effects of these nuclear burn¬ 
ing ashes on NS atmospheres and at the same time to improve 
our estimates of the NS masses and radii. 

The paper is structured as follows. In Sect. [2] we describe 
the methods we use to model the atmospheres and the emer¬ 
gent spectra. Here we present o ur new updated atmosphere code , 
originally based on the work bv ISuleimanov et ah : (1201 lb . 2012 ). 
In order to validate o ur new methods, we comp are our results to 
earlier work done by ISuleimanov et al.l ( 2012 1 and discuss the 
accuracy of the calculations. Then, in Sect. Q] a new grid of 
models is presented for various chemical compositions consist¬ 
ing of a solar mixture of H and He with highly enhanced metal- 
licities. We also determine the color correction factors for the 
metal-enriched atmospheres and show how they can be applied 
to observed X-ray burst data. At the end of this section we dis¬ 
cuss the impact these findings may have on the measurements of 
NS masses and radii and how they can act as a tool to probe the 
nuclear burning mechanisms. Finally, we summarize the main 
findings in Sect. [4] 


the modern, high-level computer programming language julia 
dBezanson et al.l [20 1 2h to allow efficient and easy scaling from 
desktop computers to supercomputing clusters. Parallelization is 
now done natively for high demand and heavy load routines like 
opacity and redistribution function calculations. 


2.1. Main equations 

In the process of modeling the atmospheres, we have to set a few 
free parameters. The first explicit input parameter in our model 
is the surface gravity 


GM , 

9=^r( 1+ z)> 


( 1 ) 


which depends on the mass M and the radius R. Here G is the 
gravitational constant and the term 1 + z= (1 - Rs/R) 12 (where 
Rs = 2GM/c 2 is the NS Schwarzschild radius and c is the speed 
of light) is the general relativistic correction. In the computations 
one must also set the stellar effective temperature T e ff of the NS. 
However, a more convenient variable is the relative luminosity 
l = L/Le dd> where L Edd is the Eddington luminosity measured at 
the surface. It is common to define F Edd as 


^Edd = 


AnGMc 

-(1 +z), 

Ke 


( 2 ) 


where the Thomson scattering opacity K e is given through the 
Thomson cross-section cr T , the gas density p and the electron 
number density N e as 


Ke = cr T — ~ 0.2(1 + X) cm 2 g *, (3) 

P 

where X is the hydrogen mass fraction. We note, however, that 
this is just a definition that we adopt to be consistent with previ¬ 
ous work. In reality, the electron cross-section is affected by the 
Klein-Nishina effect and the final form used for the Thomson 
scattering opacity K e is just an approximation where full ioniza¬ 
tion is presumed and the ratio between the number of electrons to 
protons and neutrons in the metals is assumed to be half. Defin¬ 
ing the Eddington flux and the Eddington temperature as 


^Edd - 


^Edd 

4 nR 1 


- o"sbT F hj - 


££ 

Ke ’ 


(4) 


where <x SB is the Stefan-Boltzmann constant, and noting that f = 
F/F Edd , we find a simple relation 


2. Method of atmosphere modeling 


Te« = C m T E dd. 


(5) 


The methods we use to compute the NS atmosphere mod- 
els share many import ant feat ures with earlier work s by 
ISuleimanov et al.l (1201 1 hi) and bv ISuleimanov et aD (l2012h . For 
completeness, we present the main equations and basic assump¬ 
tions in Sect. 12.11 giving emphasis to the improvements done 
for the level population and opacity calculations in Sect. 12.21 To 
solve these equations efficiently we have developed a new atmo¬ 
sphere modeling code matmos. It has its origin in the old stel¬ 
lar atmosphere modeling program ATLAS dKuruczll 1970111993 9 
which has been more recently modified to deal with high tem¬ 
perature s (llbragimov et alil2003i : ISuleimanov & Poutanenlf2005 
ISuleimanov & We rnei 12007 ) and to ta ke into account Comp¬ 
ton scattering ( Suleimanov et alJ |2012|) u sing exact relativisti c 
redistribution functions (see e.g. iPoutanen & Svenssonlll996h . 
From there on, the code was redesigned and rewritten using 


It can be used to set the effective temperature of the atmosphere 
easily via the relative luminosity. Lastly, we can freely set the 
chemical composition of our atmosphere. 

The structure of the atmosphere is described by a set of dif¬ 
ferential equations that originate from our underlying assump¬ 
tions and simplifications. The first assumption is that the atmo¬ 
sphere is in the hydrostatic equilibrium 

dP g 

-7—=g-g™d, ( 6 ) 

am 

where P g is the gas pressure, g md is the radiative acceleration and 
the column density m is found from 


dm = -pdz. 


(7) 
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where z the vertical distance. This assumption is not valid in the 
radius expansion stage of PRE bursts and, therefore, our models 
are limited to their cooling phases (and to the non-PRE bursts). 

The second important element in the theory is the radiative 
transfer equation (RTE). In our calculations we assume radiative 
equilibrium meaning that radiation field is the only means of en¬ 
ergy transport; i.e. both convection and thermal conduction are 
neglected. We can simplify the RTE by assuming a plane-parallel 
atmosphere because the scale height of the atmosphere is smaller 
(on the order of ~ 10 to ~ 10 3 cm) than the NS radius (~ 10 6 cm). 
In the plane-parallel approximation the RTE can be formulated 
with the specific intensity I(x,p) and the source function S (x, p) 
as 

dl(x,p) = i(x,n) -S(x,p), (8) 

dr(x, p) 

where p = cos 6 is the cosine of the angle between the surface 
normal and the direction of photon propagation. The dimension¬ 
less photon energy x = hv/m e c 2 is given in the units of electron 
rest mass energy m e c 2 , where h is the Planck constant, v is the 
photon frequency and m e is the electron rest mass. The optical 
depth can be related to the opacity and the column density m by 

dr(x, p) = [cr(x, p) + k(x, p)]dm, (9) 


where k is the Boltzmann constant. 

The radiation acceleration (see equation^ is then expressed 
using the RF as 


dPrad 

0rad — , — 

dm 

2n d r° r +1 2 

— — dx /rl(x,p) dp 

c dm J 0 J_i 

(14) 

-?r*j 

f [f(x, p) + k(x)] [I(x,p)-S (x, p)]pdp. 

(15) 


where the derivative with respect to m is replaced by the first 
moment of the RTE ©. The last two main equations are the 
energy balance equation (that is used to check the solution of 
RTE) 

X oo /-»+1 

J [cr(x,p) + k(x)][I(x,p)-S(x,p)]dp = 0, (16) 

and the ideal gas law (see AnncndixfAli 

Pg - NtoikT, (17) 

where N tot is the number density of all particles. In addition, par¬ 
ticle and charge conservation are implicitly assumed. 


where k(x,p ) is the true opacity due to free-free and bound- 
free transitions and cr(x,p) is the electron scattering opacity. 
Magnetic field is assumed to be negligible (the strongest mea¬ 
sured su rface mag netic field strength of a burster is < 10 10 G, 
IPanitto et al.ll201 it) in calculations of the opacity. 

In order to compute the electron scattering opacity exactly 
we need to account for induced sc attering. The exact elec tron 
scattering opacity can be written as dSuleimanov et al.ll2012l) 


i r°° r 1 

cr(x,p) — K e - I xjdxi I dp\R(x\,py,x,p) 

x Jo J-i 


1 + 


CI(x u ni) 


( 10 ) 


where R(x,p;x i,pj is the exact angle-dependent rel- 


ativistic Compton scattering redistribution function 

(RF; Aharonian & Atovan 


198lj[ Prasad et al. 

198ft 

[Nagirner & Poutanen 199A 

Poutanen & Svensson 

1996b 

that describes the probability for 

a photon with a dimensionless 


energy x propagating in a direction defined by p to be scattered 
to an energy xi and to a direction p \. The constant C is defined 
as 


C = 


1 


2 m e \m e c 2 


( 11 ) 


Now that cr(x,p) is formulated, the source function present in the 
RTE <[8]» can be written as a sum of the thermal- and the scattering 
parts 


S (x,p) 


k(x) 


cr(x,p) + k(x) 

2 


-B , + 


cr(x, p) + k(x) 


1 + 


CI(x,p) 


x .r 


r 

Jo X x J_1 


d/uiR(x,/u‘,xu/ui)I(xu/^i), ( 12 ) 


where the dimensionless Planck function B x can be written using 
the ordinary frequency dependent Planck function B v by 


dv x 3 
Bx ~ Bv dx = ~C 


1 


exp 


HtH' 


03) 


2.2. Major improvements: level population and opacity 
calculations 

Owing to the increasingly important role of metals in the plasma, 
w e have made some major im pro vements to the cod e descr ibed 
in ISuleimanov et al.1 (1201 lbl) and ISuleimanov et al.l (1201 2l) re¬ 
garding the level population and opacity calculations. We as¬ 
sume LTE and use the Saha and the Boltzmann equations in 
order to calculate the number densities of all ionization and ex¬ 
citation states, respectively. Internal partition functions are now 
build from the exact energy levels and statistical weights for all 
ions (up to Ne-like ions with a maximum net charge of X +1 °) 
and elements up to iron species with a charge of Z < 26 (ex¬ 
cluding Z = 15, 17, 19 and 21-25). Previously this was only 
done for 15 of the most abundant element^ that we have in the 
solar composition using only the ground state weights. In this 
process of building the exact partition functions, we consider 
the first 100 atomic excitation energies and statistical weights 
of states obtained from TOPbase of the Opacity Projecfl Oth¬ 
erwise the internal partition functions are built only from the 
statistical weights of the ground states like before. Neglecting 
the upper excitation levels for some heavier elements (Z > 26) 
is not, however, a great source of error due to the rapid trunca¬ 
tion of higher order terms in the internal partition function sums 
via the pressure ionization. The pressure ionization and level 
dissolution is account ed for by using the explic it occupation 
pr obability form a lism dHummer & Mihalasll 1 988l) as described 
by iHubenv et al.l (1 1994b for elements Z < 6 and otherwise th e 
method presented by the Opacity Project (ISeaton et al.lll994h . 
Previously the effects from pressure ionization was only taken 
into account with hydrogen and helium that were the main com¬ 
ponents of the atmosphere whereas in our new code it is done for 
every atomic species owing to the high concentration of heavier 
elements that can affect the general properties of the plasma. Oc¬ 
cupation probability formalism used here is an explicit (ad hoc) 
method yielding a direct correction factor that can be understood 
as a combined correction to the statistical weight (to get the real 

1 H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Fe and Ni 

2 http://cdsweb.u-strasbg.fr/topbase/home.html 
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effective weight) and to the energy level suppression due to the 
collective background electric field of the surrounding plasma. 

In addition to the exact electron scattering opacity, we take 
the free-free opacity and the bound-free transitions into account 
for all elements up to Z = 100. Again, it is important to note 
that when the fraction of metals in the composition increases, all 
of the opacity sources must be taken into account because of the 
complicated response they have for the overall opacity picture 
instead of limiting to only 15 of the most abundant elements. Es¬ 
pecially important are the elements near the iron-peak that have 
multiple bound-free edges around the spectral peak at energies 
— 1—10 keV. The bremsstrahlung opacities of all the ions are cal¬ 
culated assuming that the ion electric field is a Coulomb field of 
a charg e Ze e qual to the ion charge and using the Gaunt factors 
from lSutherland! ( 1998 ). Bound-free opacities due to photoion¬ 
ization from the ground state of all ions for elemen ts Z < 30 
is com pute d using the routine presented in IVerner & Yakovlevl 
(1 1 9951) and IVerner et al.1 d 1 9961) . For elements beyond Zn (Z = 
30) we account for the photoion ization from the hydrog en-like 
ions only using the approach by iKarzas & Latter! (119611) where 
the inner nucleus and (non-excited) electrons are replaced by a 
point-like Coulomb potential. In addition, photoionization from 
the first 20 excited levels of all hydrogen-like ions is now in¬ 
cluded using the same method. Because of the importance of the 
bound-free edges to the general properties of the emergent spec¬ 
trum it is also crucial to consider these excited levels as they, in 
many cases, smooth the otherwise sharp photoionization edges 
with additional smaller details below the ionization energy. 


2.3. Method of solution 

The calculations are performed on a logarithmic grid of 100 col¬ 
umn depth points in, ranging from 10 6 to rn max = 10 5 g cm -2 . 
Care is also taken when choosing ;« max to ensure that the effec¬ 
tive optical depth y/r v , b-f, f_f (m max )T v (m max ) is larger than unity 
at all frequencies, where r v ,b-f,f-f is the optical depth computed 
with the true opacity only (bound-free and free-free transitions, 
without scattering). This condition is necessary for satisfying the 
inner boundary condition of the radiation transfer problem (dif¬ 
fusion approximation) but it is not trivial as we need to limit 
our computational grid to the weakly coupled plasma regime 
(see Appendix [Al l. For the energy grid we use 360 logarith¬ 
mically spaced frequency points in the range 10 14 - 10 20 Hz 
(w 4x 10 4 -400 keV) for the more luminous model atmospheres 
> 0.1), and 10 14 - 10 19 Hz for l < 0.1. 

The course of the calculations is similar to that described in 
ISuleimanov et al.1 (120121) . First, as an initial guess, a gray atmo¬ 
sphere is constructed for effective temperature 7 d r (via t) from 
pre-tabulated opacity tables. Then, the radiative acceleration £/ rad 
for all depth points is ca lculated using a simple and remarkably 
accurate approximation (ISuleimanov et al.li2012l) 


(7rad 

g 


kT 

40 keV 


-l 


(18) 


The formal solution of the radiation transfer e quati on d8ji is 
obtained using the short-characteristic method (lOlson & Kunaszl 
119871) in three angles in each hemisphere. The full solution is 
then found via an acce l erated A-iteration method (see appendix 
B in ISuleimanov et ahl (1201 2l) for full description). The solution 
of the radiative transfer equation is then checked against the en¬ 
ergy balance equation (IT6l) and the surface flux condition 

r*0O 

4 71 I H x (m = 0)dx = crs B r^ ff , (19) 

Jo 


where H x is the local (astrophysical) flux at any given depth 
point m. Then the temperature corrections for every depth point 
are computed by using two different error measurements: the rel¬ 
ative error in the flux and the error in the energy balance. Here 
the relative flux error is defined 


e H (m) = 1 - 


/'CO 1 

J 0 H x (m)dx 


( 20 ) 


where Hq is the correct emergent flux corresponding to a black- 
body with a temperature 7’ e fr- In addition, the energy balance er¬ 
ror is defined as 


1 

fA (m) = - 



[cr(x, p) + k(x)][I(x, p) - S (x, p)]dp. (21) 


Corrections are then evaluated using a hybrid temperature cor¬ 
rection method consisting of three different procedures as fol¬ 
lows. In the deepest layers, the Avrett-Krook flux correction is 
used based on the relative flux error en(in). In the intermedi¬ 
ate layers, an integral A-iteration method, modified for Compton 
scattering (based on the energy balance Eq. ( IT6l) ) is used. This 
correction is done by finding the necessary temperature change 
5T in some particular depth from 


5T k = -e\(m) 

where A d (x) is the diagonal matrix element of the A-operator 
and a(x) = cn:s(x)/(k(x) + <xcs(x)) and cr cs is the Compton scat¬ 
tering opacity averaged over the rel ativistic Maxw ellian elec¬ 
tron distribution (see Eq. (A16) in lPoutanen & Svenssonl(ll996l) 
which is equivalent to Eq. (flOt if one ignores the induced scat¬ 
tering). Finally, on the uppermost layers we use the surface cor¬ 
rection based on the emerging flux error (see lKuruczlll97d for a 
detailed description). 

As a result of solving these equations iteratively, we obtain a 
self-consistent NS atmosphere model together with the emergent 
spectrum of radiation. This iteration is continued until the rela¬ 
tive flux error is less then 0.1% and the relative flux derivative er¬ 
ror is smaller then 0.01%, in most cases. For the very luminous 
models where g m dlg ~ 1 this kind of accuracy is unattainable 
and the relative flux error can be up to 1 - 2%. 


A d (x) - 1 


1 - a(x)A d (x) 


k(x } —dx 


-l 


( 22 ) 


2.4. Accuracy of computations 


that accounts for the Klein-Nishina reduction of the electron 
cross-section. The radiation pressure P rad is integrated from g rad 
yielding a good starting approximation for the pressure distribu¬ 
tions. Using these starting values, exact opacities are obtained 
for all depth and frequency points. Initial starting models built 
with the aforementioned methods serve as a good starting point 
for the iterations as the overall flux error is usually already 
smaller than 30% when compared to the actual converged model. 


Only a few metal-enriched neutron star atmospheres are present 
in the literature, so comparison and validation of our results 
is challenging. We have therefore tested some of our general 
assumption more thoroughly (see the Appendix [A}, and we 
have compared our model results of solar co mpositi on of ele¬ 
ments (SolAl) with the previous work done bv ISuleimanov et akl 
(120121) . The differences between the new- and the old models for 
the opacities obtained for solar composition with T = 10 7 K 
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(SolAl), temperature T = 10 7 K and electron pressure P = 10 13 erg 
cm - ’. The strongest bound-free edges from different chemical elements 
are indicated by arrows. 



o.l 1 10 

E (keV) 


Fig. 2. Comparison of the emergent spectra for different luminosities of 
solar composition (SolAl) computed with the current code (dash ed, col¬ 
ored lines) and the code presented in lSuleimanov et ail j2012l ) (black, 
solid lines) for log g = 14.0. Luminosities shown are L/L^m = { = 0.01 
(blue), 0.1 (red), 0.5 (green) and 1.0 (brown). The lower panel shows 
the ratio of these spectra. 


and P = 10° erg cm -3 (corresponding approximately to the 
spectral formation depth where r ~ 1) are presented in Fig.Q] 
The overall opacity picture is similar to the previous results by 
ISuleimanov et al.1 (12012l) . but as more elements are added the 
subtle details from additional bound-free edges become visible. 
Especially prominent are the additional photoionization edges 
from the excited states of hydrogen- and helium-like iron below 
the main ground state edges at 9.278 keV (Fe XXVI) and 8.828 
keV (Fe XXV). 

These small changes in the opacity profiles then mani¬ 
fest themselves as small differences in the emerging spectra 
(see Fig. 01. The results are, however, in good agreement and 
largest deviations occur in the very low-luminosity regime where 
improved level population calculations and additional opacity 
sources have the strongest effect to the emerging spectra. When 
the temperature rises and all elements become fully ionized the 
deviations due to the different edge strengths slowly vanish and 
the spectra at C = 1.0 are almost identical. The small resulting 


Table 1. Composition of the models with solar ratio of H/He and en¬ 
hanced metallicities. 


Model 

X 

Y 

Z 

A 

SolAl 

0.738 

0.249 

0.013 

1.26 

SolA20 

0.550 

0.185 

0.265 

1.65 

SolA40 

0.352 

0.118 

0.530 

2.45 


Note: X is the hydrogen mass fraction, Y is the helium mass fraction, Z 
is the metal mass fraction, and A is the mean atomic mass. 


differences are due to slightly altered electron densities. The dis¬ 
cussed differences do not have a strong impact on the previously 
derived color correction factors because the edge shapes are not 
modeled anyway by the mere diluted blackbody fit. However, 
for models with the increasing abundance of metals that act as 
sources of true opacity it becomes increasingly important to sim¬ 
ulate all the atomic species correctly. 

We also note that no bound-bound opacity is taken into 
account in the current work. Because of the partially ionized 
ions (especially H-, He- and Li-like iron ions) one indeed ex¬ 
pects some contribution from the lines to the overall opacity at 
the lower luminosities (t < 0.1 corresponding to about 7' e ir < 
20 MK). More specifically, the line blanketing would affect the 
opacity profile the most at the red side of the bound-free edges, 
rounding the sharp saw-like structures (in addition to the bound- 
free opacity of the excited states that acts similarly). In reality, 
this effect is also coupled with the rotational smearing that acts to 
smoothen out all of the high-resolution features from the spectra. 
Luckily, some computations of atmospheres with lines exists in 
this t emper ature range, namely the ones done using the TMAP 
code (IRauch et al.ll20 08) and the carb atm code (modified ver - 
sion of the code presented in this nanen lSuleimanov et al.l2014l) . 
In the end, from their results we see that the averaged continuum 
does not, however, differ much even at T e ff = 20 MK and some 
notable contributions from the lines that are strong enough to 
modify the smoothed spectra are present only at around 10 MK 
(see e.g. Fig. 8 in IRauch et al.1120081) . We also note that in prac¬ 
tice the models for X-ray bursting NS atmospheres are only used 
up to about ( ~ 0.2 but for completeness we extend our luminos¬ 
ity grid a bit farther down. 

3. New grid of models 

3. 1. General properties 

We have calculated a grid of hot NS atmosphere models using 
surface gravity values of log g — 14.0, 14.3, and 14.6, which 
cover m ost of the physically realis tic NS compactness (see e.g. 
Fig. 1 in lSuleimanov et akll201 lbi) . We considered four chemi¬ 
cal compositions: a solar ratio of H/He with metal mass fraction 
of Z = 1Z Q (SolAl) and two enhanced metallicity compositions 
with Z = 20Z o (SolA20) and Z = 40Z o (SolA40), where the 
elements from lithium to Z = 100 are enhanced by factors of 20 
and 40, respectively (see Table 1 for specific mass fractions re¬ 
garding the different models). The solar abundances were taken 
from lAsplund et al. : (120091) and for each of the compositions and 
log g values we computed models with relative luminosities l 
ranging from 0.001 to 1.06 (log gr = 14.0), 1.08 (log gr = 14.3) 
or 1.1 (log g = 14.6), depending on the chosen surface gravity. 
We have also computed a more heuristic composition consisting 
of pure iron (Fe) only, for comparison. This kind of atmosphere 
composition might not be physically realistic, but it serves as an 
extreme example of a case where the photosphere consists en- 
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tirely of heavier elements. It then acts as a proxy that sets the 
hard lower constraint for f c for a metal enriched atmosphere. 

These metal-enhanced compositions mimic the conditions in 
the photosphere of the NS in a simple way when it begins to 
cool after the rapid nuclear burning of the freshly accreted ma¬ 
terial. The key element here is the hydrogen that acts as a cata - 
lyst enabling helium burning via the op-process dSchatzll201 ll) . 
This also explains why we focus on mixed H/He bursts and not 
pure He where CNO ashes dominate. Heavy nuclei produced 
by this crp-burning then act as a seed for rp-processes produc¬ 
ing even heavier elements up to Te, far exceeding the possi¬ 
ble end-results from just 3ar-reaction or CNO-breakout ( Schat zi 
1201 ll) . The composition should still be consisting mainly of hy¬ 
drogen (and helium) but it can now be enriched with heavier 
nuclear burning ashes if there is some convection in the at¬ 
mosphere. Hydrodynamical simulations (s ee e.g. lWooslev et al.1 
120041: iFisker et all 20081: Jose et al. 1201 Oh and (semi-)analytic 


considerations (Weinberg et al, 2006) indeed show that the con¬ 


vective mixing region can approach the photosphere of the NS 
during the X-ray bursts. This can then cause an effective mixing 
of the nuclear burning ashes into the accreted metal-poor mix¬ 
ture of hydrogen and/or helium. 

The exact composition of the ashes is not yet well under¬ 
stood and many uncertainties arise from poorly known nuclear 
reaction rates, which grow into large systematic errors in the 
end-results of different nucleosynthetic simulations dParikh et al.l 
120131) . There are general trends indicating that elements around 
mass number of A ~ 20 - 30 (like Si and S) and A ~ 50 - 60 
(like Fe, Ni and Zn) have considerable overproduction factors 


when compared to solar abundances (see e.g. Koi 

te et al.l 2004 

Wooslev et al. 20041 Fisker et al. 2008; Jose et al. 

2010; Schatz 

201JJ). Because of these uncertainties, we prefer to use the sim- 


pie metal-enriched models to mimic the presence of the burn¬ 
ing ashes, rather than trying to replicate the abundances from 
the aforementioned works in our atmosphere models. The solar 
abundances of metals also peak around medium and heavy el¬ 
ements like carbon, silicon and iron, and therefore the general 
effects of these metals are likely to resemble the end-results of 
the rp-process chain (and that of moderate CNO-burning). This 
way we can find the general trends of having_metals in the mix¬ 
ture, similar to what was done bv lMaiczvna et all d2005l) where 
only iron was considered as an “average metal”. Our findings, 
however, show that it is crucial to the take a broad distribution of 
these metals into account because they not only contribute to the 
free-free continuum opacity, but also introduce bound-free edges 
that can modify the emergent spectrum considerably. Therefore, 
our models are likely to give correct general trends of having 
metals in the photosphere, which yield more electrons and con¬ 
tribute to the total free-free and bound-free opacity. 

Examples of the temperature structures and emergent spectra 
for models with log g - 14.0 and the selected chemical composi¬ 
tions are shown in Figs.[3]and[4] From the temperature structure 
of the models (see Fig. [3} one can see the general trend that more 
metal-rich atmospheres have cooler outer layersQ The tempera¬ 
ture in the upper optically thin atmospheric layers is determined 
by the balance between heating and cooling of the matter by ra¬ 
diation. The absorption opacity in the most parts of the upper 
layers is insignificant because of the low density, and the temper¬ 
ature is equal to the Compton temperature determined by heat¬ 
ing and cooling by non-coherent electron scattering only (see 


3 For the same (the metal-rich models have higher T c s because of the 

opacity dependence on X (see equations [3] and This implies that the 
effect is actually even stronger than that seen in Fig. [3] 



Fig. 3. Temperature structure for pure iron (red), solar ratio of H/He 
with 1 (black) and 40 (blue) times the solar metallicity Z e for luminosi¬ 
ties / = 0.01, 0.1 and 1.0 for log g = 14.0. 



Fig. 4. Emergent spectra for pure iron (red, dotted line), solar ratio of 
H/He with metallicity Z = 1Z Q (black solid lines) and Z = 40Z Q (blue 
dashed line) for luminosities l = 0.01, 0.1 0.5 and 1.0 computed with 
logg = 14.0. For clarity models corresponding to luminosity of l = 1.0 
are shifted along the ordinate axis by a factor of 10. 


detail discussions in lFanidus et al JlT986b iPavlov et al.lll99lh . In 
fact, the surface temperature is close to the color temperature 
of the emergent spectrum, T c = // 7’ c ft - The temperature of the 
deeper layers will be mainly determined by the balance between 

J ^OO pOO 

0 kyjydv and emission j k v B v dv, if the 
relative luminosity l is not too high. The absorption opacity in¬ 
creases toward lower photon energies and the balance occurs at 
temperatures which are much smaller than the color temperature 
of the emergent spectrum. This leads to a temperature minimum 
at some depths where the cooling and heating are dominated by 
true absorption and corresponding thermal emission. The depth 
and the width of this temperature depression depend on the ratio 
between the absorption and the scattering opacities in the atmo¬ 
sphere. Absorption is stronger in atmospheres with more abun¬ 
dant heavy elements and higher gravity (which results in a higher 
density). Indeed, the largest width of the temperature depression 
is seen for pure iron atmospheres. In the optically thick layers, 
the temperature increases according to the diffusion approxima¬ 
tion T 3 4 oc tr, where tr is the Rosseland optical depth. 

The emerging spectra are almost always harder than the cor¬ 
responding blackbody of temperature T e g (see Fig |4]). Only in 
the case of the largest surface gravity and pure iron composi¬ 
tion can the bound-free edges suppress the radiation enough, 
so that the emerging radiation becomes softer. At very high lu¬ 
minosities ( > 0.8, we confirm the previous results that the 


Article number, page 6 ofl 141 
































































J. Nattila et al.: Models of neutron star atmospheres enriched with nuclear burning ashes 


spectra become harder because of the increasing contri butio n 
from the radiation acceleration when q^aIq ~ 1 ( LondcmetdJ 
Il986t lLapidus et al .11198(3: lEbisuzakill 1987i [Pavlov et al.ll 199 ll) . 
It leads to a decrease in the matter density and the absorption- 
to-scattering opacity ratio. Therefore, the emergent photons are 
born at larger depths and at larger temperatures. At lower relative 
luminosity C (i.e. lower effective temperature), the density and 
the absorption opacity increases. This leads to the appearance 
of partially ionized species, first of all H- and He-like iron ions. 
At these lower relative luminosities, the absorption edges of the 
H-like FeXXVI (9.278 keV) and He-like FeXXV (8.828 keV) ions 
strongly affect the emerging radiation. At the coldest tempera¬ 
tures (around { ~ 0.01 - 0.1) this picture is complicated even 
further by numerous additional bound-free edges from various 
other chemical elements present in the mixture. 


3.2. Dilution and color correction factors 


The model spectra are very cl ose to a blackbody shape in the 
observed RXTE/PCA (IJahoda et al.ll2006h energy band and it is 
very common that the actual observations (especially of X-ray 
burs ts) are also fitted with a (diluted) blackbody function (see 
e.g. iGal lowav et ail 1200 8). The diluted blackbody has two pa¬ 
rameters: the color temperature and the normalization 


K = 


Rbb (km) \ = J_ / R (km) ( 1 + z) 

Dw / ic 4 V D 


'10 


(23) 


where R bb is the blackbody radius, Z>io is the distance to the 
source in units of 10 kpc and z is gravitational redshift at the 
NS surface. The theoretical evolution of f c can be related to 
the observed evolution of with flux during the cooling tail 
( Penninx et alJfl989l : Ivan Paradiis et al.lll990t ISuleimanov et al] 
1201 lallbh and the constraints on the NS mass and radius can be 
thus obtained. 

In order to make our models easily accessible for data anal- 
ysis we computed the so-called color correction curves (see e.g. 
ISuleimanov et al.l201 lbll2012l) for our new set of models. These 
can then be directly compared with observed best-fit values for 
the blackbody normalization if the em ission is coming from a 
passively cooling NS surface (see e.g. Suleimanov et~akll2012l : 
iPoutanen et alJl2014l : iKaiava et al.l 1201 41k All theoretical emer¬ 
gent spectra were fitted with a diluted blackbody 


F e ~ wB E (f c T e ff), 


(24) 


where f c is the color correction (or hardness) factor and w is 
the dilution factor. Because of the normalization of the diluted 
blackbody function we have an approximate relation f c ~ ui V 4 
between these two varia bles]^ The fitting was don e using the 
first method described in ISuleimanov et al. ( 201 lb ) in the en¬ 
ergy band (3 - 20) X (1 + z) keV corresponding to the observed 
RXTE/PCA detector energy band (see Fig. [5J- Value of the red- 
shift is calculated from the log g values by adopting NS mass of 
1.4 M q corresponding to values of z = 0.18, 0.27 and 0.42 for 
log g = 14.0 {R = 14.80 km), 14.3 (R = 10.88 km) and 14.6 
(R = 8.16 km), respectively. Evolution of the color correction 
and dilution factors are illustrated in Fig.[6]and the results of the 
fitting are presented in Table [2] for SolA20 (log g - 14.0) as an 
example. Results for other chemical compositions and gravities 
are given in the Table IB. fl 

A common feature of the color correction factors is their non¬ 
monotonic dependence on the relative luminosity t. The color 

4 wB E {f c T) dE oc wf*B(T) from where we require wf/f = 1 . 



Fig. 5. Emergent spectra (black, solid line) for different relative lu¬ 
minosities of £ = 0.01, 0.1, 0.5 and 1.0 for compositions of the solar 
ratio of H/He with Z = 2OZ 0 (upper panel) and Z = 40Z o (middle 
panel); and pure iron (bottom panel) for logg = 14.0 is shown. Best- 
fit diluted blackbody function (red, dashed line) is also shown together 
with a blackbody function where Tyt, = Tcir (blue, dotted line). Vertical 
dashed lines correspond to the observed (3 - 20) x (1 + z) keV range 
where the fitting is performed. 


correction factors have a maxima at the maximum possible (, af¬ 
ter which the correction starts to decrease as the luminosity goes 
down. As the relative luminosity decreases further, a local min¬ 
ima is finally reached around C — 0.1 - 0.6. The position of the 
minimum and its strength strongly depends on the metallicity. 
The minimum shifts toward larger (, when the heavy element 
abundance increases: for the solar abundance the minimum oc¬ 
curs at ( ~ 0.1 whereas for the pure iron atmospheres it is lo¬ 
cated at f s 0.6. This happens because the iron edges near 9 
keV suppress the high-energy radiation, forcing the most effec¬ 
tive cooling window into lower energies and thus reducing the 
color correction factor. The efficiency of this process depends on 
the occupation number of the H- and He-like iron ions. These oc¬ 
cupation numbers, on the other hand, are determined by the iron 
abundances and the ratio of the total number of iron atoms in 
different stages of ionization, which, in turn, is depended on the 
temperature and the matter density. It means that similar absorp¬ 
tion edges in the emergent spectra can arise at low temperature 
and low iron abundance as well as at a relatively high temper¬ 
ature and an increased iron abundance. This causes the shift of 
the local / c minimum to larger £, when the iron abundance is in¬ 
creased. It also explains the same track of the color corrections 
at the low luminosities for the metallicities of Z = 20Z G and 
4OZ 0 because once the edges are strong enough to suppress the 
high-energy radiation, the spectra appear rather similar. It is only 
when the peak of the effective Planck function B E (T e ff) can climb 
over this edge that the color correction starts to increase again 
with increasing L In the very low-luminosity regime (l < 0.1) 
the emerging spectra are hardened by the free-free absorption as 
first noted by iLondon et al.l d!984l 09861) . At these low photon 
energies the free-free opacity dominates over the electron scat¬ 
tering and causes additional hardening due to the strong oc £~ 3 
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Fig. 6. Color correction factors f c (solid lines) and dilution factors w 1/4 
(dashed lines) for model atmospheres consisting of pure iron (red) and 
solar mixture of H/He with Z = Z G (SolAl, black), Z = 20Z G (SolA20, 
magenta) and Z = 40Z o (SolA40, blue) against the relative luminosity l 
for logt/ = 14.0 (upper panel), log g = 14.3 (middle panel) and log g = 
14.6 (bottom panel). 


energy dependency. The effect is then further enhanced by the 
decreasing effective temperature T e n leading to an even harder 
spectrum. 

Another notable feature is the apparent super-Eddington lu- 
mi nosities (when the relativ e luminosity is f > 1) as first noted 
bv lSuleimanov et £ 01 ( 120121 ) . This happens because of the Klein- 
Nishina reduction of the electron cross-sections. Because we 
compute our models all the way up to g md /g ~ 1, formally super- 
Eddington luminosities are needed to counter the fact that the 
Eddington luminosity is defined for the Thomson cross-section. 
As Klein-Nishina reduction depends on the electron temperature 
it also makes the upper limiting luminosity dependent on log g 
as the models with stronger surface gravity tend to be hotter. We 
also note that in our definition of the Eddington luminosity (and 
hence also in the Eddington temperatures tabulated) we use the 
approximation K e ~ 0 . 2(1 + X ) cm 2 g 1 for the Thomson scatter- 


Table 2. Color correction and dilution factors from the blackbody fits 
to the spectra of SolA20 atmosphere models at log g = 14.0. 


Z = 20Z o log g = 14.0 T Edd = 1.75 keV 
R= 14.80 km z = 0.18 


{ 

gmd/g 

T eff (keV) 

fc 

w 

wf? 

0.001 

0.017 

0.311 

1.805 

0.002 

0.023 

0.003 

0.026 

0.409 

1.508 

0.023 

0.121 

0.010 

0.043 

0.553 

1.286 

0.185 

0.505 

0.030 

0.067 

0.728 

1.229 

0.371 

0.847 

0.050 

0.087 

0.827 

1.233 

0.401 

0.927 

0.070 

0.106 

0.900 

1.235 

0.416 

0.968 

0.100 

0.139 

0.984 

1.232 

0.438 

1.008 

0.150 

0.190 

1.089 

1.214 

0.484 

1.049 

0.200 

0.240 

1.170 

1.192 

0.529 

1.067 

0.300 

0.332 

1.295 

1.177 

0.546 

1.046 

0.400 

0.414 

1.391 

1.201 

0.496 

1.033 

0.500 

0.526 

1.471 

1.238 

0.415 

0.975 

0.550 

0.568 

1.507 

1.258 

0.386 

0.968 

0.600 

0.610 

1.540 

1.280 

0.359 

0.964 

0.650 

0.651 

1.571 

1.302 

0.335 

0.962 

0.700 

0.692 

1.600 

1.326 

0.312 

0.962 

0.750 

0.734 

1.628 

1.351 

0.290 

0.966 

0.800 

0.775 

1.655 

1.378 

0.270 

0.972 

0.850 

0.816 

1.680 

1.404 

0.251 

0.975 

0.900 

0.857 

1.704 

1.435 

0.232 

0.982 

0.940 

0.890 

1.723 

1.464 

0.216 

0.990 

0.960 

0.907 

1.732 

1.478 

0.208 

0.991 

0.980 

0.921 

1.741 

1.494 

0.199 

0.993 

1.000 

0.939 

1.750 

1.519 

0.188 

1.000 

1.020 

0.953 

1.758 

1.551 

0.173 

0.998 

1.040 

0.972 

1.767 

1.599 

0.154 

1.008 

1.060 

0.989 

1.775 

1.634 

0.140 

1.012 


Notes. Results for other chemical compositions and surface 
gravities are listed in the Appendix iBl 


ing opacity. This approximation assumes full ionization and that 
the ratio of number of electrons to protons and neutrons in metals 
is half, which, in reality, is not the case. In practice, these correc¬ 
tions are, however, rather small and easy to apply, so for clarity 
and convenience we use the standard definition of the Eddington 
luminosity defined by equations (0 and 0. 

The new models also show significant dependency on the 
surface gravity because of the metals. When increasing the log g 
value the density of the photosphere is also enhanced. This in¬ 
crease in density, in turn, leads to smaller iron ionization frac¬ 
tions and, therefore, to larger bound-free edges as the number 
density of absorbing elements increase. The luminosity where 
the large iron edges are finally seen to vanish (due to almost full 
ionization) can be traced to £ = 0.5 for SolA20 and C = 0.6 for 
SolA40 with log g = 14.0 whereas the same limiting luminosity 
appear around C - 0.8 for SolA20 and f — 0.9 for SolA40 when 
the surface gravity is increased to log g - 14.3. The disappear¬ 
ance of these edges then lead to a break in the color correction 
curves when photons are able to escape from the blue side of the 
photoionization edge. This change results in a sudden increase in 
the hardness of the emerging spectra and hence increase in the 
/ c . For the log g = 14.6 this break is shifted to near the Edding¬ 
ton luminosity and coincides with the color correction increas¬ 
ing when the luminosity becomes close to the Eddington limit. 
These two effects in turn make the color correction evolution 
for models with log g = 14.6 very sharp at ( > 0.8. Moreover, 
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the color correction evolutions t — f c for both metal-enriched 
compositions are very similar when log g - 14.6. The similar¬ 
ity originates because the increased temperature with SolA40 
model (that will reduce the number of absorbing nuclei) is coun¬ 
tered by the larger initial number of iron nuclei. As a result 
the model spectra of SolA20 and SolA40 became self-similar 
as the SolA20 model, in contrast, has a lower temperature (i.e. 
more non-ionized nuclei that can absorb) but fewer heavy ele¬ 
ment nuclei (like iron) capable of participating in the absorption 
processes. The ( — f c curves for pure iron models do not show 
any significant increase because iron is not fully ionized even at 
l ~ 1 . 

From our color correction fits we can limit the variability 
of / c to be between about 1.3 - 1.8 near the Eddington limit and 
1.0-1.4 around t ~ 0.5. At lower luminosities of around t ~ 0.1 
the scatter is even smaller (/ c ~ 1.2-1.4) but there is a possibility 
that at this late stage the spectra of the cooling NS surface are 
already likely to be contaminated by additional heating from the 
in-falling gas from the accretion flow. It should also be noted that 
the lower limit set by the pure iron composition is most likely 
not realistic as at least some hydrogen (and other elements too 
from the ashes of the previous bursts) is likely to be present in 
the atmosphere. Still, it acts as a heuristic limiting case, with a 
maximum concentration of metals. It is also interesting to note 
that even small amounts of metals are able to modify the color 
correction factor substantially. 

It is also clear that the spectra of NS with pure iron envelopes 
and even just strongly enhanced heavy element abundances can¬ 
not be fitted well with smooth blackbody functions. Therefore, 
our computations of color factors for their model spectra are 
rather arbitrary. In fact, when comparing to observations we 
would expect to find significant residuals near iron absorption 
edges if the atmosphere indeed has large (heavy) metal content. 
However, we have to also keep in mind that various reasons, 
such as the fast rotation of the NS, additional spectral smear¬ 
ing because of the rapidly rotating (spreading) boundary layer or 
additional radiation from the optically thin recombination con- 
tinua of iron ions arising in a surrounding envelope can reduce 
the observed strength of the iron absorption edges. 

In addition to the direct observations of the iron edges, the 
connection between the color correction factor and the black- 
body normalization (given by equation [23l> opens up a possibil¬ 
ity of observing the presence of these metals indirectly: even a 
small increase in their abundance would end up decreasing the 
color correction factor toward unity. Our results show that this 
mechanism is enough to explain the small burst-to-burst vari¬ 
ations observed in the coolin g tracks between separate bursts 
from the same source (see e.g. Giiver et akll2012tlPoutanen et al.1 
l2014t> . On the other hand, abnormally large deviations from the 
typical cooling track for some individual burst might indicate an 
especially large ash content and/or effective mixing in the pho¬ 
tosphere. In addition, these outlier bursts are expected to show 
strongest residuals (i.e. poor x 2 values) when the photosphere 
cools giving yet another testable prediction of the models. These 
connections to the observations remain to be tested and will be 
pursued in future publications. 

4. Summary 

We have presented a new NS atmosph ere modeling code which 
is a n extension of our previous works (iSuleimanov et al.ll2011 hi 
120121) . The code is especially designed to model the emerging 
spectra from X-ray bursting NSs with metal-rich atmospheres 
and can be used for NS mass and radius estimates. The compu¬ 


tational scheme and the code are validated producing results in a 
good agreement with our previous work. 

We have assumed an ideal plasma in LTE. The model atmo¬ 
spheres are considered to be in hydrostatic and radiative equi¬ 
librium. We use the standard plane-parallel atmosphere approx¬ 
imation. No magnetic field nor the stellar rotation are taken into 
account. For Compton scattering, the exact angle-dependent rel¬ 
ativistic redistribution kernel is used. As a new addition, the 
bound-free and free-free opacities are now taken into account for 
almost all (up to Z = 100) elements. We also use exactly con¬ 
structed internal partition functions for calculations of the ion¬ 
ization fractions and use the occupation formalism to take the 
pressure ionization effects into account in comparison to previ¬ 
ous work. 

To study the effects arising from the metals we have com¬ 
puted a new set of hot NS atmospheres with enhanced metallici- 
ties. Atmospheric structures and emergent spectra are calculated 
for four different chemical compositions (pure iron and solar hy¬ 
drogen/helium mix with various abundances of heavy elements 
of Z = 1, 20 and 40 Z Q ) and for three different surface gravities 
of log g - 14.0, 14.3 and 14.6. The relative luminosities range 
from C = 0.001 to 1.06-1.10 (depending on log g). 

All computed emergent spectra differ from the blackbody 
spectra. The spectral shape for the relatively luminous models 
is, however, well reproduced by a diluted blackbody that is used 
to fit all of our models in the observed RXTE/PCA energy band 
(3 - 20 keV)@ The spectra for the less luminous models have 
significant iron absorption edges, and their approximation by a 
diluted blackbody becomes worse. 

We found that the emergent spectra are strongly dependent 
on the fraction of (heavy) metals in the composition. We con¬ 
strain the color correction factor to be between f c ~ 1.3 - 1.8 
near ( » 1 and / c ~ 1.0 - 1.4 around t ~ 0.5, depending on 
the metallicity. These limits show that by varying the fraction of 
the burning ashes in the composition, the color correction can 
change quite considerably. It also shows that - in addition to the 
uncertain detection of the photoionization edges in the emergent 
spectra - the color correction evolution can be used as a new tool 
to probe the nuclear burning (and the mixing via convection) at 
the NS surface due to the strong metal dependency. 
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Optical depth r R 

Fig. A.l. Plasma coupling parameter T as a function of Rosseland 
mean optical depth tr. Gray area marks the regime of (strong) coupling 
where our original assumption of ideal gas law slowly breaks. Presented 
are coupling parameters for atmosphere compositions of pure iron (red 
lines) and solar mixture of H and He with metallicity Z = 1Z Q (SolAl, 
black lines) and Z = 40Z o (SolA40, blue lines) computed for hydrogen 
(solid lines), carbon (dashed) and iron (dot-dashed) ions. 


Appendix A: Validity of ideal plasma assumption 

Compared to previous atmosphere models computed with hydro¬ 
gen or helium compositions our work differs substantially due to 
high number densities of metals producing large electron num¬ 
ber densities. In the very deep layers, the dense plasma might be¬ 
come strongly coupled invalidating our assumption of ideal gas 
law (O- To quantify the depth where the coupling starts to be 
important, we can compare the ion specific Coulomb energies 
Ec and local thermal energies Erh to get the so-called plasma 
coupling parameter, which is defined as 


E c _ Z 2 e 2 
Ej h r s kT ’ 


(A.l) 


/ i \ 1/3 

where r s = (^L-1 is the Wigner-Seitz radius, i.e. the typical 
inter-particle distance, Z is the charge of the ion and n e is the 
local electron density (see e.g. lCheniri984l for basic plasma pa¬ 
rameters). For an ideal plasma (with equation of state that of the 
ideal gas), the thermal (kinetic) energy of ions is larger than the 
interaction potentials. This corresponds to the situation T <K 1. 
In the opposite regime, pair distribution of ions becomes strongly 
localized and the plasma starts to crystallize into an (ideal) solid 
state. Between these two states there exists a mixed phase of 
fluid-like plasma around T ~ 1. The area where our approach is 
valid is limited to the case of the weakly coupled plasma, T <s 1. 

The T-values are presented in Fig. lA.ll for different compo¬ 
sitions and ions. The regime of intermediate coupling F ~ 1 is 
reached only for heavier ions in the electron-rich dense layers 
at Rosseland optical depths > 10 4 . Thus this cannot affect the 
emergent spectra. However, this can cause problems if the max¬ 
imum column depth m max is chosen to be too large so that the 
regime F > 1 is reached within our computational domain. In 
these very deep layers the inter-particle effects may become sub¬ 
stantial, leading to a numerically challenging instability between 
the strong (and unphysical) recombination via the Saha equation 
and forced ionization because of the occupation probability for¬ 
malism used. This is a known problem because Saha’s model 
assumes ideal gas conditions (i.e. every effect of the plasma 
is asserted on its ions and atoms internal structures) and uses 
only the ground state c onfigurations in the ionization balance 
equations ([Rogers et al.lll996 ;). It is also worth mentioning that 
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the pressure ionization method used here is a more of a phe¬ 
nomenological correction scheme as its corrections to non-ideal 
effe cts are obtained by minimi zing a free energy of predefined 
gas dHummer & Mihalasll988lf . More fundamental and physical 
approach would be the (activity) expansion of the grand canoni¬ 
cal partition function of the plasma where pressure ionization is 
a consequence of the theory (see e.g. iRogers et all! 19961) . Such 
a fundamental treatment of plasma is, however, extremely chal¬ 
lenging and out of the scope of this paper. 

Appendix B: Color correction tables 

Table B.l. Color correction and dilution factors from the blackbody fits. 
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0.158 

0.949 

0.300 

0.312 

1.539 

1.147 

0.576 

0.998 

z 

= 40 Z o 

log g = 14 

.0 r Edd = 1.81 

keV 

0.400 

0.390 

1.654 

1.182 

0.493 

0.964 


R ■■ 

= 14.80 km 

z = 0.18 


0.500 

0.471 

1.748 

1.226 

0.420 

0.951 

0.001 

0.026 

0.322 

1.503 

0.006 

0.032 

0.550 

0.510 

1.791 

1.248 

0.391 

0.949 

0.003 

0.041 

0.424 

1.335 

0.048 

0.151 

0.600 

0.548 

1.830 

1.269 

0.366 

0.950 

0.010 

0.064 

0.573 

1.228 

0.239 

0.544 

0.650 

0.589 

1.867 

1.292 

0.342 

0.952 

0.030 

0.093 

0.754 

1.213 

0.394 

0.851 

0.700 

0.633 

1.902 

1.315 

0.320 

0.956 

0.050 

0.113 

0.857 

1.222 

0.415 

0.926 


Article number, page 12 ofl 141 














J. Nattila et al.: Models of neutron star atmospheres enriched with nuclear burning ashes 


( 

0rad / 9 

T eS (keV) 

fc 

w 

wj* 

t 

(7rad / 9 

Tes (keV) 

fc 

w 

wf c 4 

0.070 

0.131 

0.932 

1.226 

0.429 

0.967 

1.060 

0.966 

2.185 

1.535 

0.177 

0.984 

0.100 

0.163 

1.019 

1.223 

0.451 

1.008 

1.080 

0.971 

2.196 

1.623 

0.141 

0.982 

0.150 

0.214 

1.128 

1.206 

0.498 

1.053 

Z 

= 40Z o 

log g = 14 

•6 r Edd - 2.56 

keV 

0.200 

0.263 

1.212 

1.180 

0.557 

1.079 


R 

= 8.16 km 

z = 0.42 


0.300 

0.354 

1.341 

1.141 

0.640 

1.084 

0.001 

0.022 

0.455 

1.859 

0.011 

0.125 

0.400 

0.428 

1.441 

1.126 

0.654 

1.053 

0.003 

0.032 

0.599 

1.474 

0.093 

0.440 

0.500 

0.505 

1.524 

1.133 

0.614 

1.012 

0.010 

0.048 

0.809 

1.272 

0.327 

0.855 

0.550 

0.543 

1.561 

1.143 

0.585 

0.999 

0.030 

0.074 

1.065 

1.212 

0.482 

1.040 

0.600 

0.599 

1.595 

1.168 

0.531 

0.988 

0.050 

0.103 

1.210 

1.165 

0.599 

1.104 

0.650 

0.671 

1.627 

1.208 

0.451 

0.960 

0.070 

0.129 

1.317 

1.123 

0.714 

1.134 

0.700 

0.726 

1.657 

1.241 

0.402 

0.954 

0.100 

0.164 

1.439 

1.080 

0.835 

1.135 

0.750 

0.763 

1.686 

1.266 

0.371 

0.953 

0.150 

0.215 

1.593 

1.054 

0.876 

1.081 

0.800 

0.800 

1.714 

1.293 

0.341 

0.954 

0.200 

0.252 

1.712 

1.066 

0.833 

1.074 

0.850 

0.837 

1.740 

1.323 

0.312 

0.957 

0.300 

0.372 

1.894 

1.111 

0.621 

0.946 

0.900 

0.873 

1.765 

1.357 

0.284 

0.963 

0.400 

0.455 

2.036 

1.162 

0.509 

0.928 

0.940 

0.900 

1.784 

1.387 

0.261 

0.967 

0.500 

0.535 

2.152 

1.211 

0.432 

0.929 

0.960 

0.913 

1.794 

1.404 

0.249 

0.970 

0.550 

0.575 

2.204 

1.233 

0.403 

0.932 

0.980 

0.931 

1.803 

1.426 

0.237 

0.979 

0.600 

0.614 

2.253 

1.254 

0.378 

0.937 

1.000 

0.943 

1.812 

1.447 

0.224 

0.980 

0.650 

0.653 

2.298 

1.275 

0.357 

0.942 

1.020 

0.960 

1.821 

1.476 

0.209 

0.994 

0.700 

0.693 

2.341 

1.294 

0.338 

0.949 

1.040 

0.974 

1.830 

1.513 

0.191 

0.999 

0.750 

0.732 

2.382 

1.313 

0.322 

0.955 

1.060 

0.989 

1.839 

1.559 

0.170 

1.004 

0.800 

0.770 

2.421 

1.332 

0.306 

0.963 

z 

= 40Z o 

logs' = 14.: 

3 r Edd = 2.15 

keV 

0.850 

0.809 

2.458 

1.350 

0.292 

0.968 


R 

= 10.88 km 

z = 0.27 


0.900 

0.846 

2.493 

1.369 

0.277 

0.973 

0.001 

0.024 

0.383 

1.621 

0.010 

0.071 

0.940 

0.877 

2.520 

1.388 

0.264 

0.979 

0.003 

0.036 

0.504 

1.362 

0.091 

0.313 

0.960 

0.892 

2.534 

1.398 

0.257 

0.982 

0.010 

0.055 

0.681 

1.221 

0.341 

0.759 

0.980 

0.905 

2.547 

1.409 

0.249 

0.982 

0.030 

0.080 

0.896 

1.231 

0.412 

0.945 

1.000 

0.924 

2.560 

1.427 

0.239 

0.993 

0.050 

0.101 

1.018 

1.223 

0.451 

1.009 

1.020 

0.938 

2.572 

1.444 

0.228 

0.992 

0.070 

0.126 

1.108 

1.208 

0.493 

1.049 

1.040 

0.950 

2.585 

1.469 

0.213 

0.991 

0.100 

0.160 

1.211 

1.179 

0.563 

1.087 

1.060 

0.958 

2.597 

1.504 

0.193 

0.986 

0.150 

0.209 

1.340 

1.131 

0.679 

1.112 

1.080 

0.971 

2.609 

1.570 

0.163 

0.988 

0.200 

0.253 

1.440 

1.103 

0.745 

1.102 

1.100 

0.981 

2.621 

1.665 

0.129 

0.991 

0.300 

0.327 

1.594 

1.082 

0.762 

1.046 


Fe log g = 14.0 

r Edd = 

1.95 keV 

0.400 

0.395 

1.713 

1.089 

0.707 

0.993 


R ■■ 

= 14.80 km 

z = 0.18 


0.500 

0.465 

1.811 

1.114 

0.622 

0.958 

0.001 

0.022 

0.347 

1.998 

0.001 

0.014 

0.550 

0.500 

1.855 

1.131 

0.579 

0.947 

0.003 

0.047 

0.457 

2.061 

0.007 

0.130 

0.600 

0.533 

1.895 

1.147 

0.542 

0.940 

0.010 

0.081 

0.618 

1.836 

0.048 

0.548 

0.650 

0.564 

1.934 

1.164 

0.510 

0.935 

0.030 

0.182 

0.813 

1.587 

0.137 

0.866 

0.700 

0.595 

1.970 

1.181 

0.480 

0.933 

0.050 

0.290 

0.924 

1.486 

0.198 

0.969 

0.750 

0.624 

2.004 

1.198 

0.453 

0.933 

0.070 

0.370 

1.005 

1.425 

0.248 

1.024 

0.800 

0.654 

2.037 

1.215 

0.429 

0.934 

0.100 

0.435 

1.099 

1.365 

0.309 

1.072 

0.850 

0.683 

2.068 

1.232 

0.406 

0.936 

0.150 

0.521 

1.216 

1.304 

0.385 

1.114 

0.900 

0.724 

2.098 

1.256 

0.375 

0.935 

0.200 

0.588 

1.307 

1.252 

0.464 

1.139 

0.940 

0.851 

2.121 

1.294 

0.318 

0.892 

0.300 

0.692 

1.446 

1.148 

0.667 

1.160 

0.960 

0.860 

2.132 

1.333 

0.293 

0.925 

0.400 

0.770 

1.554 

1.091 

0.814 

1.152 

0.980 

0.897 

2.143 

1.398 

0.256 

0.976 

0.500 

0.833 

1.643 

1.060 

0.880 

1.112 

1.000 

0.924 

2.154 

1.431 

0.233 

0.977 

0.550 

0.857 

1.683 

1.055 

0.875 

1.083 

1.020 

0.946 

2.164 

1.458 

0.218 

0.987 

0.600 

0.880 

1.720 

1.055 

0.850 

1.054 

1.040 

0.964 

2.175 

1.490 

0.202 

0.996 

0.650 

0.898 

1.754 

1.062 

0.809 

1.028 
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l 

gmdlg 

T eS (keV) 

fc 

w 

wJt 

0.700 

0.915 

1.787 

1.072 

0.759 

1.003 

0.750 

0.929 

1.818 

1.088 

0.701 

0.983 

0.800 

0.946 

1.848 

1.108 

0.641 

0.967 

0.850 

0.954 

1.876 

1.134 

0.578 

0.956 

0.900 

0.970 

1.903 

1.165 

0.514 

0.948 

0.940 

0.971 

1.924 

1.197 

0.461 

0.947 

0.960 

0.982 

1.934 

1.214 

0.436 

0.947 

0.980 

0.984 

1.944 

1.233 

0.410 

0.948 

1.000 

0.987 

1.954 

1.254 

0.385 

0.950 

1.020 

0.988 

1.964 

1.273 

0.363 

0.953 

1.040 

0.989 

1.973 

1.298 

0.338 

0.959 

1.060 

0.991 

1.983 

1.329 

0.309 

0.964 


Fe 

log g — 14.3 

?Edd = 

2.32 keV 


R 

' = 10.88 km 

z = 0.27 


0.001 

0.025 

0.413 

2.117 

0.004 

0.072 

0.003 

0.045 

0.543 

2.029 

0.021 

0.359 

0.010 

0.074 

0.734 

1.727 

0.085 

0.753 

0.030 

0.199 

0.967 

1.472 

0.213 

1.001 

0.050 

0.307 

1.098 

1.369 

0.307 

1.080 

0.070 

0.365 

1.194 

1.309 

0.382 

1.120 

0.100 

0.428 

1.306 

1.253 

0.468 

1.154 

0.150 

0.510 

1.445 

1.170 

0.633 

1.185 

0.200 

0.576 

1.553 

1.098 

0.820 

1.192 

0.300 

0.678 

1.719 

1.031 

1.024 

1.158 

0.400 

0.757 

1.847 

1.006 

1.060 

1.084 

0.500 

0.818 

1.953 

1.006 

0.991 

1.014 

0.550 

0.843 

2.000 

1.013 

0.936 

0.986 

0.600 

0.865 

2.044 

1.024 

0.875 

0.964 

0.650 

0.885 

2.085 

1.040 

0.809 

0.946 

0.700 

0.902 

2.124 

1.059 

0.742 

0.933 

0.750 

0.918 

2.161 

1.082 

0.675 

0.925 

0.800 

0.932 

2.196 

1.108 

0.611 

0.920 

0.850 

0.945 

2.230 

1.138 

0.549 

0.920 

0.900 

0.956 

2.262 

1.171 

0.490 

0.924 

0.940 

0.964 

2.287 

1.201 

0.446 

0.929 

0.960 

0.968 

2.299 

1.217 

0.425 

0.933 

0.980 

0.972 

2.311 

1.234 

0.404 

0.937 

1.000 

0.976 

2.322 

1.252 

0.384 

0.942 

1.020 

0.979 

2.334 

1.270 

0.364 

0.947 

1.040 

0.982 

2.345 

1.290 

0.344 

0.953 

1.060 

0.985 

2.356 

1.312 

0.324 

0.960 

1.080 

0.988 

2.367 

1.336 

0.304 

0.968 


Fe 

log 0 = 14.6 

?Edd = 

2.76 keV 


= 8.16 km 

z = 0.42 


0.001 

0.025 

0.491 

2.080 

0.013 

0.238 

0.003 

0.041 

0.646 

1.905 

0.045 

0.598 

0.010 

0.068 

0.873 

1.573 

0.152 

0.929 

0.030 

0.214 

1.149 

1.319 

0.373 

1.129 

0.050 

0.306 

1.305 

1.223 

0.530 

1.187 

0.070 

0.356 

1.420 

1.168 

0.654 

1.215 


e 

{7rad / 9 

r eff (keV) 

fc 

w 

wf c 4 

0.100 

0.416 

1.552 

1.090 

0.875 

1.236 

0.150 

0.496 

1.718 

1.010 

1.175 

1.222 

0.200 

0.560 

1.846 

0.975 

1.297 

1.171 

0.300 

0.663 

2.043 

0.949 

1.296 

1.052 

0.400 

0.740 

2.195 

0.956 

1.157 

0.965 

0.500 

0.799 

2.321 

0.982 

0.983 

0.914 

0.550 

0.824 

2.377 

1.001 

0.896 

0.899 

0.600 

0.846 

2.429 

1.023 

0.814 

0.891 

0.650 

0.866 

2.478 

1.047 

0.737 

0.887 

0.700 

0.885 

2.525 

1.073 

0.668 

0.886 

0.750 

0.902 

2.569 

1.101 

0.605 

0.890 

0.800 

0.917 

2.610 

1.130 

0.550 

0.896 

0.850 

0.931 

2.650 

1.159 

0.500 

0.904 

0.900 

0.944 

2.688 

1.190 

0.456 

0.914 

0.940 

0.954 

2.718 

1.214 

0.424 

0.923 

0.960 

0.959 

2.732 

1.228 

0.409 

0.928 

0.980 

0.963 

2.746 

1.241 

0.394 

0.934 

1.000 

0.968 

2.760 

1.254 

0.379 

0.939 

1.020 

0.972 

2.774 

1.268 

0.365 

0.945 

1.040 

0.976 

2.787 

1.283 

0.351 

0.951 

1.060 

0.980 

2.801 

1.299 

0.336 

0.958 

1.080 

0.984 

2.814 

1.317 

0.321 

0.966 

1.100 

0.987 

2.827 

1.336 

0.306 

0.974 
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